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A new kind of overset grid, named Yin- Yang grid, for spherical geometry is proposed. The Yin- 
Yang grid is composed of two identical component grids that are combined in a complemental way 
. . . to cover a spherical surface with partial overlap on their boundaries. Each component grid is a 

, low latitude part of the latitude-longitude grid. Therefore the grid spacing is quasi-uniform and 

■ the metric tensors are simple and analytically known. One can directly apply mathematical and 
' numerical resources that have been written in the spherical polar coordinates or latitude-longitude 

grid. The complemental combination of the two identical component grids enables us to make 
. efficient and concise programs. Simulation codes for geodynamo and mantle convection simulations 

using finite difference scheme based on the Yin- Yang grid are developed and tested. The Yin- Yang 
\ grid is suitable for massively parallel computers. 
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I. INTRODUCTION 

. . 

' ' Since the Earth is composed of spherical layers, computer simulations of the Earth's interior, such as geodynamo and 
^ \ mantle convection simulations, need efficient spatial discretization schemes in spherical she ll g eometry. The spectral 
method fisj has been the major tool in the geodynamo simulation; all six codes jfll Ull li^ l4(l l47|| in the benchmark 
^ ' test in Christensen et al. flOj and other codes [e.g., JQ, 2S] use the spherical harmonics expansion method in the 
^ t horizontal space. However, the importance of non-spectral (or point-based) approaches in the dynamo simulation 
' ^ ■ is now increasingly recognized to simulate more realistic geodynamo regime with smaller Ekman numbers • The 
pursuit of point-based approaches started earlier in the mantle convection simulations, because the mantle's intense 
spatial variation of viscosity and the phase transitions makes the spectral approach not fit to the problem. Although 
Qh, the spectral method for the mantle convection prospered in 1980s and 90s (161, [13, 13^ .4^] , the finite element method 

■ is rapidly growing in this field [Hp , l3^ liM l50l | . There are also a couple of codes that uses the finite element method in 
' the geodynamo simulation Ia l8ll | . The finite difference or finite volume method is applied for the mantle convection 

^ . by Hernlund and Tackley [l9| , Iwase [2H , Ratcliff et al. ,38.J . The finite difference method has been used for the core 
rn ' convection and the geodynamo simulation by the authors from 1990s 0, |23, 113, 13> which the 

, latitude- longitude grids in the spherical polar coordinates is used with radius r (r^ < r < Tq), colatitude 9 (0 < 9 < tt), 
' and longitude 4> {0 < (f) < 27r). Since the finite difference method enables us to make highly optimized programs for 
massively parallel computers, especially massively parallel vector supercomputers like the Earth Simulator |l7j |. we 
further exploit the possibility of the finite difference method for simulations in spherical shell geometry by improving 
[ the base grid system. 

, It is known that the latitude-longitude grid has two numerical problems; the coordinate singularity and the grid 
' convergence near the poles. Since the coordinate singularity is not a real singularity (the pole is not singular point of 
, physical functions), one can solve the basic equations on the poles by applying the I'Hospital's rule on the pole grids 
J>^" [e-g-jHil- The computational cost for this pole grid solver is negligible. 

] The problem of the grid convergence is more serious. In order to relax the severe restriction on the time step, 
Oh' one has to apply a filter so that the grid spacing on the sphere becomes effectively quasi- uniform. The amount of 
^ ' information abandoned by the filter is estimated by the number of grid points that are effectively present and that 
actually present in the computational space; suppose one has a latitude-longitude grid of a spherical surface of unit 
radius with inter mesh angles A in both colatitude (9) and longitude (</>). The azimuthal grid spacing, which is A in 
the equator, converges in higher latitudes. When a filter enables an effectively quasi-uniform grid with spacing A on 
the sphere, the number of effective grid points is estimated by {J^ sin 9 d9 J^'^ d4>)/A^ = 47r/A^. While the number 
of actual grid points in the computational space is given by {J^ d9 x J^^ (i(/))/A^ — 27r^/A^. Therefore sizeable ratio 
of information, (27r^ — Att)/2tt'^ ^ 36% of the latitude- longitude grid, is abandoned in vain by the filtering at each 
simulation step. In addition to this computational inefficiency, the filter has non-negligible computational costs. In 
our geodynamo simulation code using latitude-longitude grid, in which a Fast Fourier Transform (FFT)-based filtering 
procedure is applied, the filter routine can take more than 30% of the total execution time. 
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Note that the above problem of the grid redundancy in the latitude-longitude grid comes only from the region of 
high latitudes. The remaining part of the latitude-longitude grid — the low latitude region-has rather desirable feature 
for numerical simulations; it is an orthogonal grid, it has simple metric tensors, and it has quasi-uniform grid spacings. 
This observation leads us to the idea of a new spherical grid proposed in this paper. 

Since there is no grid mesh that is orthogonal all over the spherical surface and, at the same time, free of coordinate 
singularity or grid convergence, we decompose the spherical surface into subregions. The decomposition, or dissection, 
enables us to cover each subregion by a grid system that is individually orthogonal and singularity-free. This divide- 
and-rule approach has been used with good success in the computational aerodynamics that incorporates complex 
geometry of aircraft's body with wings/stores/blades. 

The dissection of the computational domain generates internal border or boundary between the subregions. There 
are two different approaches to handle the internal boundaries. One is the patched grid method [s^l and the other is 
the overset grid method |^. In the patched grid approach, the subdomains contact one another without any overlap 
on their borders. In the overset grid method, on the other hand, the subdomains partially overlap one another on 
their borders. The overset grid is also called as overlaid grid, or composite overlapping grid, or Chimera grid [43 |. 
The validity and importance of the overset approach in the aerodynamical calculations was pointed out by Steger 
|43| . Since then this method is widely used in this field. It is now one of the most important grid techniques in the 
computational aerodynamics; for ex amp le, whole aircraft with wing and store |3^ . tiltrotor aircraft |33j . Boeing 747 
1^1 Hoi 1 Space Shuttle 5], helicopter and others. 

In the computational geosciences, the idea of the overset grid approach appeared rather early. Phillips proposed 
a kind of composite grid in 1950's to solve partial differential equations on a hemisphere, in which the high latitude 
region of the latitude-longitude grid is "capped" by another grid system that is constructed by a stereographic 
projection to a plane on the north pole |3l l35ll3a |. After a long intermission, the overset grid method seems to attract 
growing interest in geoscience these days. The "cubed sphere" is an overset grid that covers a spherical surface 
with six component grids that correspond to six faces of a sphere. The "cubed sphere" is recently applied to the 
mantle convection simulation In the atmospheric research, other kind of spherical overset grid is used in a global 
circulation model [T^ . in which the spherical surface is covered by two component grids — improved stereographic 
projection grids — in northern and souther hemispheres that overlap in the equator. A successful test of 100-day 
integration of global circulation is demonstrated with this overset grid. 

The overset grid proposed in this paper is named "Yin- Yang grid" after the symbol for yin and yang of Chinese 
philosophy of complementarity. The Yin- Yang grid is composed of two identical and complemental component grids. 
Compared with other spherical overset grids, the Yin- Yang grid is simple in its geometry and metric tensors. A 
remarkable feature of this overset grid is that the two identical component grids are combined in a complemental way 
with a special symmetry. 



II. BASIC YIN- YANG GRID 



The Yin- Yang grid in its most basic shape is shown in Fig. ^ It has two component grids that are geometrically 
identical (exactly the same shape and size); see Fig. QJa). We call the two component grids "Yin grid" (or n-grid) 
and "Yang grid" (or e-grid). They are combined to cover a spherical surface with partial overlap on their borders as 
shown in Fig. ^b). Each component grid is in fact a part of the latitude-longitude grid: A component grid, say Yin 
grid, is defined in the spherical polar coordinates by 

(7r/4 -6<9< 37r/4 + S) n (-37r/4 - S < (/) < Sir/A + S), (1) 

where J is a small buffer, which is proportional to grid spacing, required for minimum overlap in the overset method- 
ology (see Fig. ^b)). In the limit of infinitesimal grid {6 — > 0), the area of the above part of the sphere with unit 

radius is given by J^J^^ sin 9d9 J'^^^J^^d(j> = 3tt/\/2 ~ 2.127r, i.e., roughly a half of the whole spherical surface (27r). 

Another component grid, Yang grid, is defined by the same rule of eq. ^ but in different spherical coordinates that 
is perpendicular to the original one; see the green- and blue-colored spherical mesh in Fig. |21 The axis of the Yang 
grid's coordinates (blue mesh in Fig.|21l, is located in a equator of the Yin grid's coordinates (green mesh in Fig.l^J. 
The relation between Yin coordinates and Yang coordinates is denoted in the Cartesian coordinates by 

(a;^y^z^) = (-x^z^y"), (2) 

where (a;", i/", z") is Yin's Cartesian coordinates and (x*^, y'^, z^") is Yang's. In a matrix form, 

(3) 
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where 

M = I 1 I . (4) 
Note that 

Af-i = A/, (5) 

which indicates that the transformations between Yin and Yang coordinates are symmetric. This is a reflex of the 
complemental relation between Yin and Yang. 
In the spherical coordinates, eq. (jSJ reads 

= r", (6) 

sinrcos(/)*= = -sinr cos(/)", (7) 

sin sin 0"^ = cose*", (8) 

cosO" = sin r sin 0", (9) 

where (r", 6'", (/)"), and {r^^O^^cIf) are the coordinates of Yin and Yang, respectively. The idea of two perpendicular 
spherical coordinates is used in the global ocean simulation to avoid the grid convergence in the Arctic, however, 
the second spherical coordinates is used in a sort of auxiliary way for the main (usual) spherical polar coordinates in 
their method. On the other hand, we make the best use of the symmetry between two coordinates. 

For spatial discretization, we define mesh point at j-th colatitude and /c-th longitude cf/p. on Yin grid (for I = n) 
and on Yang grid (for I — e) as 

- Vin + J (j = 0, A^e - 1), (10) 
= 0niin + ^^'^' {k^O,N^-l), (11) 

with 

= (^max- Vin)/(^e-l), (12) 

A0 = (0max-0niin)/(^^- 1)' (13) 

where the grid distribution ranges from 9^:^^ = 7r/4 — 5 to 6'niax — 37r/4 + (5 in colatitude, and from (/^j-qJ^ = — 37r/4 — (5 
to (/>max = 37r/4 + 5 in longitude. We set KO — Acj) = 25 in Fig.Q] as an example. 

An important feature of the Yin- Yang grid as a spherical overset grid is that the two component grids are identical 
and their geometrical positions are complemental. This enables us to make concise programs: Suppose a grid point 
{6",(j)^j^) on Yin grid's horizontal border at index position {j,k) (e.g., j = 1). Its value should be determined by an 
interpolation from its neighbor points, or stencils, of Yang grid with interpolation coefficients that are determined by 
relative position of (6*",^^) in the stencils. Note that exactly the same interpolation coefficients and relative stencils 
are used to set the value of corresponding grid point (0|, c/)^) at {j, k) of Yang's border, since the geometrical relations 
between Yin grid and Yang grid are symmetric. In other words, we can make use of one interpolation routine for two 
times (for Yin grid and for Yang grid) to set the horizontal boundary conditions. Note also that the metric tensors at 
a bulk grid point at (j, k) of Yin grid is a function of its position (0", cf)^) in Yin's coordinates, and the metric tensors 
at corresponding point (0|, <f>'j.) in Yang grid are exactly the same. Therefore we can call one subroutine of fluid solver 
and others for two times for Yin grid and Yang grid. 

Another advantage of the Yin- Yang grid resides in the fact that the component grid is nothing but the (part of ) 
latitude-longitude grid. We can directly deal with the equations to be solved with the vector form in the usual spherical 
polar coordinates, {w^, t'e, i'^}- The analytical form of metric tensors are familiar in the spherical coordinates. We 
can directly code the basic equations in the program as they are formulated in the spherical coordinates. We can 
make use of various resources of mathematical formulas, program libraries, and tools that have been developed in the 
spherical polar coordinates. 

To conclude this section, we point out that the construction of three-dimensional Yin- Yang grid for spherical shell 
geometry is straightforward, by piling up the basic (two-dimensional) Yin- Yang grids in radial direction. See Fig. |31 



III. VECTOR TRANSFORMATION FORMULA BETWEEN YIN AND YANG GRIDS 



Following the general overset methodology [e.g., 8], interpolations are applied on the boundary of each component 
grid to set the boundary values, or internal boundary condition. When one deals with scalar variables, the interpolation 
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is simple. For vector fields, a care is needed for vector components, since expressions of a vector in the Yin's spherical 
coordinates, {VrjVgjV"^}, and in the Yang's coordinates, {vf. , Vg , v^} , are difi'erent. 

Because the Yin- Yang transformation denoted by eq. lO is a rotation about the origin (r = 0), the radial component 
of the vector is invariant (uj? = v^), and horizontal components are mapped by local rotation transforms, as shown in 
Fig. 21 where the rotation angle is a function of latitude and longitude; 

,n \ 
r \ 

(14) 

To find the expression of ■0, we consider unit vectors in 9 and (p directions on the Yin and Yang coordinates. From 
Fig. ^ we see 

cosV = (15) 
sinV- = (16) 

where r and 0^ are unit vectors in 6 and (p directions in the component grid £, with £ = n for Yin grid, and £ — e for 
Yang grid. The unit vectors {x^, y^, z^} in the Cartesian coordinates are related to {9^, 0^} by 

= -sin</>'^F +cos0^y^ (17) 

0" = -sin</-"5:" +cos<?!." J/" 

= sin<^!)"^'' + cos0"z^ (18) 

r = cos 9" cos (p^x" + 003 9" sin -81119" z''. (19) 

Substituting eqs. (|17|l and (|18|) into (|15|l . we get 

cos ip — — sin (j)'^ sin </>" . (20) 
Substituting eqs. and into we get 

sin ip — cos 6**^ cos (j)'^ sin 0" — sin 9"^ cos 

^ {cos9'' (sin 6*'= cos 0*= ) (sin 61" sin 0" ) - sin^ 6*^ (sin 9" cos (f)'' ) } 



sm W'- sm 

COS( 



sin 6'" 

cos (/)" 

sin 6*'= ■ 



(21) 



Here we have used eqs. (01 -(jU. 

From eqs. (|^ . l(T5|l . (fT?^ and ifT^ . we obtain the transformation formula of the vector components (wr, V0,v^) 

by 

/ V'^\ ( V'] \ 

U| = P , (22) 

\-%) Vil 

with the transformation matrix 

10 

-sin(/>'=sin(/)" - cos sin ^ | . (23) 
cos 0" / sin 9"^ — sin 0'^ sin 0" 

Since Yin and Yang coordinates are symmetric, the inverse transformation from Yang into Yin is given by the 
interchange of the suffixes: 

10 

P^^ = I -sin0"sin0'= -cos0Vsin6'" | . (24) 
cos 4>" / sin 6*" — sin sin < 
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Note also that 

P2 ^ 1, (25) 

which indicates the complemental relation between Yin and Yang coordinates. 

When we see the component grid of the basic Yin- Yang grid shown in Fig. ^ in the Mercator projection, it is a 
rectangle; the four corners intrude most into the other component grid (see Fig. H^b)). Even if the grid mesh is 
taken to be infinitesimal, i.e., A9 = Acj) — > and 5 — > 0, the overlapping area has still non-zero ratio of about 6%; 
{3/V2tt - 27r)/27r ~ 0.061. This overlapped area can be minimized by modifying the component grid's shape from 
the rectangle. It is obvious that a Yin- Yang grid with minimum overlap region can be constructed by a division, or 
dissection, with a closed curve on a sphere that cuts the sphere into two identical parts. There are infinite number of 
such dissections of a sphere. Fig. |S1 shows two examples among them. When we cut along the curve that is colored 
with red and blue in Fig.|3fa) or (b), we get two separated parts of the spherical surface that are identical. Although, 
it is not apparent that the two parts separated by the blue-red curve in each panel of Fig. [S] are identical from this 
figure, the corresponding three-dimensional view (Fig.EJ would show more convincingly. The cutoff curve of Fig.|^a) 
reminds us a baseball, while the cutoff curve of Fig. Etb) resembles a cube. 

Based on these spherical dissections, we can construct spherical overset grids with two identical component grids 
that has minimum overlapping area; Fig. [3 shows a Yin- Yang grid that corresponds to the baseball type dissection 
of a sphere (panels labeled (a) in Figs. Inland IHJl. Fig. |S|is for the cube type dissection (panels (b) of Figs. |31and|nj). 
When minimizing the computational cost is strongly required, the Yin- Yang grid of the baseball type (Fig. [71) or cube 
type (Fig. ISJ would be worth trying. 

However, the non-rectangle geometries of the component grid of Fig. [7|or Fig. |S1 imply that special cares should be 
taken to mask some grid points. The number of the mask is the same for both the Yin- Yang grids of Figs. [7| and |H1 
since the non-masked area of a component grid is just a half of the spherical surface {2tt) in the limit of the negligibly 
small overlap area. 

IV. SUMMARY 

For numerical simulations of the Earth's interior, we have developed a new spherical grid based on the overset grid 
methodology. Our motivation is to devise an spherical grid system that is suitable for finite difference scheme on 
massively parallel vector supercomputers. The spherical overset grid proposed in this paper, named Yin- Yang grid, 
is composed of two component grids. They have the same shape and size and combined to cover a spherical surface 
with partial overlap on their borders. Each component grid is nothing but low latitude region of the usual latitude- 
longitude grid; it is 90° about the equator and 270° in the longitude. Therefore the grid spacing is quasi-uniform and 
the metric tensors are simple and analytically known. One can directly apply mathematical and numerical resources 
that have been written in the spherical polar coordinates or latitude-longitude grid system. Since the two component 
grids are identical and combined in a complemental way, various routines for solvers and interpolation can be recycled 
for two times for each component grid at every simulation time step. 

We have developed finite difference codes of the mantle convection and dynamo simulation using the basic Yin- Yang 
grid for spherical shell geometry (see Figs.n^^iidlSl). We have confirmed that the Yin- Yang grid is successfully applied 
to both cases. The mantle convection code is newly developed from scratch. Details of the code and simulation 
results are reported in other paper |48j |: we solved the time development of thermal convection motion in a spherical 
shell of a Boussinesq fluid with infinite Prandtl number for uniform and variable viscosity cases. We have performed 
standard benchmark tests of the mantle convection js^, and confirmed that the results of our Yin- Yang mantle 
convection code successfully reproduced previously published results. The numerical values of Nusselt number and 
the mean velocity coincides with other benchmark values within a few percent or even better We have also 

applied the Yin- Yang grid to the geodynamo simulation code. The magnetohydrodynamic (MHD) equations with 
finite viscosity, thermal diffusivity, and electrical conductivity are solved. The Yin- Yang geodynamo code has been 
converted from our previous geodynamo code which was based on the latitude-longitude grid. We found that the 
code conversion was straightforward and rather easy since the base grid is common. We could reproduce our previous 
(latitude-longitude grid based) results of geodynamo simulation by our newly developed Yin- Yang geodynamo code 
with shorter calculation time. The details of the code will be reported in other paper. 

The Yin- Yang grid is suitable for parallel programming. Since the number of the component grid is two, we are 
naturally lead to make parallel programs with domain decomposition of even number: We first decompose whole 
computational region into two — Yin component and Yang component — then apply further domain decomposition in 
each component. 

Finally, we point out another possible spherical overset grid that has an odd number of component grids. Fig. |51 
shows a spherical overset grid that consists of three identical component grids. In this case, the component grid is 
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defined as a part (about 1/3) of tlie splierical surface by (7r/4 < 9 < 37r/4) fl {—n/2 < 4> < 7r/2). This grid could be 
effective when the processor number is multiple of three. 
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FIG. 1: Basic Yin- Yang grid, (a) It is a spherical overset grid composed of two identical component grids, named Yin and 
Yang, (b) The Yin grid and Yang grid are combined to cover a spherical surface with partial overlap. 




FIG. 3: Three-dimensional Yin- Yang grid for spherical shell geometry. This is constructed by piling up the basic Yin- Yang 
grid shown in Fig. Q in the radial direction. 




FIG. 4: Unit horizontal vectors of Yin and Yang coordinates. They are mapped one another by a local rotation transform with 
angle Tp. 
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FIG. 5: Curves that divide a spherical surface into two identical areas. If one cuts along the blue-red curve of the panel (a) or 
(b), the spherical surface is divided into two parts (denoted by Yin and Yang in the pictures) that are exactly the same shape 
and size. The blue part of the curve and red part of the curve are in the complemental relation; The blue curve of Yin is red 
curve of Yang, and vice versa. 




FIG. 6; Curves that divide a spherical surface into two identical areas. These are corresponding three-dimensional views of the 
red-blue curves in (a) and (b) of Fig. |S] 




FIG. 7; A Yin- Yang grid with minimum overlap that has the baseball-hke border curve between Yin and Yang grids. Corre- 
sponding spherical dissection is Figs. |^a) audita). 




FIG. 8: Another Yin- Yang grid with minimum overlap. The border curve between Yin and Yang grids is cube-like, 
sponding spherical dissection is Figs. I^b) andl^b). 



Corre- 




FIG. 9: Another possible spherical overset grid that is composed of three identical component grids. 



